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Abstract. - Based on Sirovich's two-fluid kinetic theory and a dodecagonal discrete velocity 
model, a two-dimensional 61-velocity finite-difference lattice Boltzmann method for the com- 
plete Navier-Stokes equations of binary fluids is formulated. Previous constraints, in most 
existing lattice Boltzmann methods, on the studied systems, like isothermal and nearly in- 
compressible, are released within the present method. This method is designed to simulate 
compressible and thermal binary fluid mixtures. The validity of the proposed method is ver- 
ified by investigating (i) the Couette flow and (ii) the uniform relaxation process of the two 
components. 



Introduction. - Lattice Boltzmann Method (LBM) has become a viable and promising 
numerical scheme for simulating fluid flows. There are several options to discretize the Boltz- 
mann equation: (i) Standard LBM (SLBM) [1]; (ii) Finite-Difference LBM (FDLBM) [1-3]; 
(hi) Finite- Volume LBM [1,4]; (iv) Finite- Element LBM [1,5]; etc. These kinds of schemes 
are expected to be complementary in the LBM studies. 

Even though various LBMs for multicomponent fluids [6-19] have been proposed and 
developed , (i) most existing methods belong to the SLBM [6-16], and/or based on the single- 
fluid theory [7-14,16,17,20]; (ii) in Ref. [6] a SLBM based on Sirovich's two-fluid kinetic theory 
[21] is proposed; (hi) nearly all the studies are focused on isothermal and nearly incompressible 
systems. In a recent study [22], Sirovich's kinetic theory is clarified and corresponding two- 
fluid FDLBMs for Euler equations and isothermal Navier-Stokes equations are presented. In 
this letter we propose a two-fluid FDLBM for the complete Navier-Stokes equations, including 
the energy equation. 

Formulation and verification of the FDLBM. - The formulation of a FDLBM consists of 
three steps: (i) select or design an appropriate discrete velocity model (DVM), (ii) formulate 
the discrete local equilibrium distribution function, (hi) choose a finite-difference scheme. 
The continuous Boltzmann equation has infinite velocities, so the rotational invariance is 
automatically satisfied. Recovering rotational invariant macroscopic equations from a discrete- 
finitc-vclocity microscopic dynamics imposes constraints on the isotropy of DVM and the 
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finite-difference scheme used. In this Letter, the proposed FDLBM is based on the following 
DVM, 



v = 0, v fci ; = v k 
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where k indicates the fc-th group of particle velocities and i indicates the direction of the 
particle speed. It is easy find that (i) its odd rank tensors are zero, and (ii) its initial four 
even rank tensors satisfy 
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where a, (3, ■ ■ ■ indicate x or y component and 
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It is clear that this DVM is isotropic up to, at least, its 9th rank tensor. 

We consider a binary mixture with two components, A and B, where the masses and 
temperatures of the two components are not significantly different. The interparticle collisions 
can be divided into two kinds: collisions within the same species (self-collision) and collisions 
among different species (cross-collision) [21]. Based on the DVM the 2-dimensional BGK 
[23] kinetic equation for species A reads [22], 
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j A (o) anc i jAB(o) are |- ne corresponding Maxwellian distribution functions. n A , u A , T A are 
the local density, hydrodynamic velocity and temperature of species A. u AB , T AB are the 
hydrodynamic velocity and temperature of the mixture after equilibration process. a A is the 
acceleration of species A due to the effective external field. 
For species A, we have 
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where P A (e^ t ) is the local pressure (internal mean kinetic energy). For species B, we have 
similar relations. For the mixture, we have 
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where p A = n A m A , n = n A + n B and p = p A + p B . Three sets of hydrodynamic quantities 
(for the two components A, B and for the mixture) are involved, but only two sets of them 
are independent. So this is a two-fluid model. Without lossing generality, we focus on hy- 
drodynamics of the two individual species. By expanding the local equilibrium distribution 
function f AB (°) around f A (°> to the first order in flow velocity and temperature, the BGK 
model (|(jl9l) becomes 
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where p A = p A p B /(T AB p), p A = k B n A n B /(r AB n), M = n A p A p B / \2t ab np) . 

Now, we go to the second step: formulate f^°^- The continuous Maxwellian f A ^ pos- 
sesses an infinite sequence of moment properties. The Chapman-Enskog analysis [24] shows 
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Recall that u A (u B ) is a small quantity. By using Eq. I|15|l . P A — n A kBT A , and neglecting 
the second and higher order terms in u A , Eq. Ijltjfl shows that the diffusion velocity, — u A , 
is related to the gradients of n A and T A . 

,A(0) 



The first three requirements on f ki (U> are referred to Eq. I|1U|) with f ki replaced by fuf°\ 



and the remaining five are 
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The requirement equation l|23[) contains the fifth order of the flow velocity u A . 
sufficient to expand f k }°^ in polynomial up to the fifth order of u A : 
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So it is 
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The truncated equilibrium distribution function f k ^ (|24l) contains the fifth rank tensor of 
the particle velocity w A and the requirement l|20l) contains its third rank tensor. Thus, a 
DVM being isotropic up to its 8th rank tensors is enough to recover the physical isotropy of 
the continuous Boltzmann equations to the Navier-Stokes level. So DVM is an appropri- 
ate choice. To calculate the discrete f ki , one first needs calculate the factor F A . F A is 
determined by the eight requirements on f ki and the isotropic properties of the DVM 
Following the same procedure as described in [22], we obtain 
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Once a zero speed, v A = 0, and other five nonzero ones, v A (k = 1, 2, 3, 4, 5) are chosen, F A 
(k = 0, 1, 2, 3, 4, 5) will be fixed. 

We come to the third step: finite-difference implementation of the discrete kinetic method. 
There are several choices [17] available. One possibility is shown below, 
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where the second superscripts n, n + 1 indicate the consecutive two iteration steps, At the 
time step; the spatial derivatives are calculated as 

dfti (n) = / (3/ fc 1;f } - + f#%)/{2Aa) if vi a > 

^ I (3/,1f - + /^i" ) 2 )/(-2Aa) if < Q < ' 

where a = x,y, the third subscripts J — 2, J — 1, J, I + 1, 7 + 2 indicate consecutive mesh 
nodes in the a direction. 

The validity of the formulated FDLBM is verified through two test examples. (The Boltz- 
mann constant kg — 1.) The first one is the isothermal and incompressible Couette flow with 
a single component. In this case, A = B. The initial state of the fluid is static. The distance 
between the two walls is D. At time t = they start to move at velocities U, —U, respectively. 
The horizontal velocity profiles of species A or B along a vertical line agree with the following 
analytical solution, 

u = ,V- E(-D J+1 ^ exp(-^t) sin(^), (29) 

where 7 = 2U/D is the imposed shear rate, j is an integer, the two walls locate at y = ±D/2. 
(For example, see Fig. 1.) 

The second one is the uniform relaxation process, which is an ideal process to indicate the 
equilibration behavior of the mixture [22] . By neglecting the force terms and terms in spatial 
derivatives, the Navier-Stokes equations l|15l) - (|17l) give 

= 0= (30) 
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The flow velocities of the two components equilibrate exponentially with time. (For example, 
see Fig. 2(a).) The equilibration of flow velocities also affects that of the temperatures. When 
the flow velocity difference is zero, the temperatures equilibrate exponentially with time. (For 
example, see Fig. 2(b).) The simulation results agree well with Eqs. I|31|l and (|32|l . 

Conclusions and remarks. - The Chapman-Enskog analysis shows what properties the 
discrete Maxwellian distribution function should follow. Those requirements tell the 

lowest order of the flow velocity u A in the Taylor expansion of f^°^- The highest rank of 
tensors of the particle velocity v A in the requirements on the truncated f^°^ determines 
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y/D t t 



Fig. 1 Fig. 2 

Fig. 1 - Horizontal velocity profiles along a vertical line for the two species, A and B, at time t = 8. 
The symbols are for simulation results. The solid line corresponds to the theoretical result, Eq. 1291 . 
Parameters used in the two-fluid FDLBM are m A = m B = 1, T = 1, n A = n B = 1, 7 = 0.001, 
t aa _ t bb _ t ab = t ba = Q2 p arameters use d in gq. are 77 = r) A = 0.1, p = p A = 1. 



Fig. 2 - Uniform relaxation processes, (a) Equilibration of velocities; (b) Equilibration of tem- 
peratures. The symbols are for simulation results. The solid lines possess the theoretical slopes. 
Common parameters for the simulations in (a) and (b) are n A = 10, n B = 1, m A — 1, m B = 10, 

T AA = r BB = 1) T AB = 1Q) T BA = ± j n ^ ^ j^j^ conditions are U A (°) = _ M f (°) = „ .3, 

uf 0) = u B (0) = 0, and T A(0) = 1.3, T s(0) = 0.7. The slope of the solid line in (a) is -11/20, which 
is consistent with Eq. |gTJ. In (b) the initial conditions are u A(0) = u s(0) = 0, and T A(0) = 1.3, 
yB(o) _ 0.7. The slope of the solid line in (b) is —10.1/11, which is consistent with the first term 
of right-hand side of Eq. 11321 . The second superscript "(0)" denotes the corresponding initial value. 
This figure shows an example where the particle masses of the two species are significantly different. 



the needed isotropy of the DVM. The incorporation of the force terms makes no additional 
requirement on the isotropy of the DVM. The present approach works for binary neutral 
fluid mixtures. One possibility to introduce interfacial tension is to modify the pressure 
tensors [13], which is implemented by changing the force terms [3]. The specific force terms 
or pressure tensors, which are out of the scope of this Letter, depend on the system under 
consideration, but can be resolved under the same two-dimensional 61-velocity model (D2V61). 
For binary fluids with disparate-mass components, say m A <C m B , only if the total masses and 
temperatures of the two species are not significantly different, does Sirovich's kinetic theory 
works [22], so do the corresponding FDLBMs. (See Fig. 2 for an example.) When the masses 
and/or the temperatures of the two components are greatly different, the two-fluid kinetic 
theory should be modified. In those cases, the Navier-Stokes equations and the FDLBMs 
are not symmetric about the two components, but the FDLBMs can still be resolved under 
the D2V61 model. The formulation procedure is straightforward. In practical simulations, 
numerical errors from the finite-difference schemes should be quantified. 



* * * 
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